!> \file GFS_suite_interstitial_4.F90
!!  Contains code to calculate tendencies of tracers due to convective transport, updates tracers after convective transport, and updates cloud condensation nuclei.

  module GFS_suite_interstitial_4

  contains

!> \section arg_table_GFS_suite_interstitial_4_run Argument Table
!! \htmlinclude GFS_suite_interstitial_4_run.html
!!
    subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, tracers_total, ntrac, ntcw, ntiw, ntclamt, &
      ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, ntccn, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson,  &
      imp_physics_nssl, nssl_invertccn, nssl_ccn_on,                                                  &
      imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, dtidx, dtend,&
      index_of_process_conv_trans, gq0, clw, prsl, save_tcp, con_rd, con_eps, nssl_cccn, nwfa, spechum, ldiag3d,     &
      qdiag3d, save_lnc, save_inc, ntk, ntke, otsptflag, errmsg, errflg)

      use machine,               only: kind_phys
      use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber

      implicit none

      ! interface variables

      logical, intent(in)     :: otsptflag(:)! on/off switch for tracer transport by updraft and
      integer,              intent(in   )                   :: im, levs, tracers_total, ntrac, ntcw, ntiw, ntclamt, ntrw, &
        ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, ntccn, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson,    &
        imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imp_physics_nssl

      logical,                                  intent(in) :: ltaerosol, convert_dry_rho
      logical,                                  intent(in) :: nssl_ccn_on, nssl_invertccn

      real(kind=kind_phys), intent(in   )                   :: con_pi, dtf
      real(kind=kind_phys), intent(in   ), dimension(:,:)   :: save_qc
      ! save_qi is not allocated for Zhao-Carr MP
      real(kind=kind_phys), intent(in   ), dimension(:,:)   :: save_qi, save_lnc, save_inc

      ! dtend and dtidx are only allocated if ldiag3d
      logical, intent(in)                                   :: ldiag3d, qdiag3d
      real(kind=kind_phys), dimension(:,:,:), intent(inout) :: dtend
      integer,              dimension(:,:),   intent(in)    :: dtidx
      integer,                                intent(in)    :: index_of_process_conv_trans,ntk,ntke

      real(kind=kind_phys), dimension(:,:,:), intent(inout) :: gq0
      real(kind=kind_phys), dimension(:,:,:), intent(inout) :: clw
      real(kind=kind_phys), dimension(:,:),   intent(in) :: prsl
      real(kind=kind_phys),                   intent(in) :: con_rd, con_eps, nssl_cccn
      real(kind=kind_phys), dimension(:,:),   intent(in) :: nwfa, save_tcp
      real(kind=kind_phys), dimension(:,:),   intent(in) :: spechum

      character(len=*),     intent(  out)                   :: errmsg
      integer,              intent(  out)                   :: errflg

      ! local variables
      real(kind=kind_phys), parameter :: zero = 0.0_kind_phys, one = 1.0_kind_phys
      integer :: i,k,n,tracers,idtend
      real(kind=kind_phys) :: liqm, icem, xccn, xcwmas, xccw, xcimas, qccn

      real(kind=kind_phys) :: rho, orho
      real(kind=kind_phys), dimension(im,levs) :: qv_mp !< kg kg-1 (dry mixing ratio)
      real(kind=kind_phys), dimension(im,levs) :: qc_mp !< kg kg-1 (dry mixing ratio)
      real(kind=kind_phys), dimension(im,levs) :: qi_mp !< kg kg-1 (dry mixing ratio)
      real(kind=kind_phys), dimension(im,levs) :: nc_mp !< kg-1 (dry mixing ratio)
      real(kind=kind_phys), dimension(im,levs) :: ni_mp !< kg-1 (dry mixing ratio)

      ! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      ! This code was previously in GFS_SCNV_generic_post, but it really belongs
      ! here, because it fixes the convective transportable_tracers mess for Zhao-Carr
      ! and GFDL MP from GFS_suite_interstitial_3. This whole code around clw(:,:,2)
      ! being set to -999 for Zhao-Carr MP (which doesn't have cloud ice) and GFDL-MP
      ! (which does have cloud ice, but for some reason it was decided to code it up
      ! in the same way as for Zhao-Carr, nowadays unnecessary and confusing) needs
      ! to be cleaned up. The convection schemes doing something different internally
      ! based on clw(i,k,2) being -999.0 or not is not a good idea.
      do k=1,levs
        do i=1,im
          if (clw(i,k,2) <= -999.0) clw(i,k,2) = 0.0
        enddo
      enddo

      if(ldiag3d) then
         if(ntk>0 .and. ntk<=size(clw,3)) then
            idtend=dtidx(100+ntke,index_of_process_conv_trans)
            if(idtend>=1) then
               dtend(:,:,idtend) = dtend(:,:,idtend) + clw(:,:,ntk)-gq0(:,:,ntk)
            endif
         endif
         if(ntcw>0) then
            if (imp_physics == imp_physics_zhao_carr     .or. &
                 imp_physics == imp_physics_zhao_carr_pdf .or. &
                 imp_physics == imp_physics_gfdl) then
               idtend=dtidx(100+ntcw,index_of_process_conv_trans)
               if(idtend>=1) then
                  dtend(:,:,idtend) = dtend(:,:,idtend) + clw(:,:,1)+clw(:,:,2) - gq0(:,:,ntcw)
               endif
            else if(ntiw>0) then
               idtend=dtidx(100+ntiw,index_of_process_conv_trans)
               if(idtend>=1) then
                  dtend(:,:,idtend) = dtend(:,:,idtend) + clw(:,:,1)-gq0(:,:,ntiw)
               endif
               idtend=dtidx(100+ntcw,index_of_process_conv_trans)
               if(idtend>=1) then
                  dtend(:,:,idtend) = dtend(:,:,idtend) + clw(:,:,2)-gq0(:,:,ntcw)
               endif
            else
               idtend=dtidx(100+ntcw,index_of_process_conv_trans)
               if(idtend>=1) then
                  dtend(:,:,idtend) = dtend(:,:,idtend) + clw(:,:,1)+clw(:,:,2) - gq0(:,:,ntcw)
               endif
            endif
         endif
      endif

!  --- update the tracers due to deep & shallow cumulus convective transport
!           (except for suspended water and ice)

      if (tracers_total > 0) then
        tracers = 2
        do n=2,ntrac
!         if ( n /= ntcw .and. n /= ntiw .and. n /= ntclamt) then
!          if ( n /= ntcw  .and. n /= ntiw  .and. n /= ntclamt .and. &
!               n /= ntrw  .and. n /= ntsw  .and. n /= ntrnc   .and. &
!               n /= ntsnc .and. n /= ntgl  .and. n /= ntgnc  &
!               .and. &
!             n /= nthl  .and. n /= nthnc .and. n /= ntgv    .and. &
!             n /= nthv .and. n /= ntccn  &
!                                                               ) then
           IF ( otsptflag(n) ) THEN                                                    
              tracers = tracers + 1
            if(n/=ntk .and. n/=ntlnc .and. n/=ntinc .and. n /= ntcw .and. n /= ntiw) then
               idtend=dtidx(100+n,index_of_process_conv_trans)
               if(idtend>=1) then
                  dtend(:,:,idtend) = dtend(:,:,idtend) + clw(:,:,tracers)-gq0(:,:,n)
               endif
            endif
            do k=1,levs
              do i=1,im
                gq0(i,k,n) = clw(i,k,tracers)
              enddo
            enddo
          endif
        enddo
      endif

      if (ntcw > 0) then

!  for microphysics
        if (imp_physics == imp_physics_zhao_carr     .or. &
            imp_physics == imp_physics_zhao_carr_pdf .or. &
            imp_physics == imp_physics_gfdl) then
           gq0(1:im,:,ntcw) = clw(1:im,:,1) + clw(1:im,:,2)

        elseif (ntiw > 0) then
          do k=1,levs
            do i=1,im
              gq0(i,k,ntiw) = clw(i,k,1)                     ! ice
              gq0(i,k,ntcw) = clw(i,k,2)                     ! water
            enddo
          enddo

          if ( imp_physics == imp_physics_nssl ) then
              liqm =  con_pi/6.*1.e3*(18.e-6)**3  ! 4./3.*con_pi*1.e-12
              icem =  con_pi/6.*1.e3*(120.e-6)**3 ! 4./3.*con_pi*3.2768*1.e-14*890.
              qccn = nssl_cccn/1.225 !1.225 is a reference air density and should match what is used in module_mp_nssl_2mom.F90 (rho00)
              do k=1,levs
                do i=1,im
                   ! check number of available ccn
                   IF ( nssl_ccn_on ) THEN
                     IF ( nssl_invertccn ) THEN
                       xccn = qccn - gq0(i,k,ntccn)
                     ELSE
                       xccn = gq0(i,k,ntccn)
                     ENDIF
                   ELSE
                     xccn = Max(0.0, qccn - gq0(i,k,ntlnc))
                   ENDIF
                   
                   IF ( gq0(i,k,ntlnc) > 0.0 .and. save_qc(i,k) > 0.0 ) THEN
                      xcwmas = Max( liqm, clw(i,k,2)/gq0(i,k,ntlnc) )
                   ELSE
                      xcwmas = liqm
                   ENDIF

                   IF ( gq0(i,k,ntinc) > 0.0 .and. save_qi(i,k) > 0.0 ) THEN
                      xcimas = Max( liqm, clw(i,k,1)/gq0(i,k,ntinc) )
                   ELSE
                      xcimas = icem
                   ENDIF
                   
                  IF ( xccn > 0.0 ) THEN
                  xccw = Min( xccn, max(0.0, (clw(i,k,2)-save_qc(i,k))) / xcwmas )
                  gq0(i,k,ntlnc) = gq0(i,k,ntlnc) + xccw 
                  IF ( nssl_ccn_on ) THEN
                     IF ( nssl_invertccn ) THEN
                       ! ccn are activated CCN, so add
                       gq0(i,k,ntccn) = gq0(i,k,ntccn) + xccw
                     ELSE
                       ! ccn are unactivated CCN, so subtract
                       gq0(i,k,ntccn) = gq0(i,k,ntccn) - xccw
                     ENDIF
                  ENDIF
                  ENDIF

                  gq0(i,k,ntinc) = gq0(i,k,ntinc)  &
                           +  max(0.0, (clw(i,k,1)-save_qi(i,k))) / xcimas
                enddo
              enddo
          endif

          if (imp_physics == imp_physics_thompson .and. (ntlnc>0 .or. ntinc>0)) then
            if_convert_dry_rho: if (convert_dry_rho) then
              do k=1,levs
                do i=1,im
                  !> - Convert specific humidity to dry mixing ratio
                  qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k))
                  !> - Density of air in kg m-3 and inverse density
                  rho = con_eps*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+con_eps))
                  orho = one/rho
                  if (ntlnc>0) then
                    !> - Convert moist mixing ratio to dry mixing ratio
                    qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k))
                    !> - Convert number concentration from moist to dry
                    nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k))
                    nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho)
                    !> - Convert number concentrations from dry to moist
                    gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k))
                  endif
                  if (ntinc>0) then
                    !> - Convert moist mixing ratio to dry mixing ratio
                    qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k))
                    !> - Convert number concentration from moist to dry
                    ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) 
                    ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho)
                    !> - Convert number concentrations from dry to moist
                    gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k))
                  endif
                enddo
              enddo
            else
              do k=1,levs
                do i=1,im
                  !> - Density of air in kg m-3 and inverse density
                  rho = con_eps*prsl(i,k) / (con_rd*save_tcp(i,k)*(spechum(i,k)+con_eps))
                  orho = one/rho
                  if (ntlnc>0) then
                    !> - Update cloud water mixing ratio
                    qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k))
                    !> - Update cloud water number concentration
                    gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho)
                  endif
                  if (ntinc>0) then
                    !> - Update cloud ice mixing ratio
                    qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k))
                    !> - Update cloud ice number concentration
                    gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho)
                  endif
                enddo
              enddo
            end if if_convert_dry_rho
            if(ldiag3d .and. qdiag3d) then
              idtend = dtidx(100+ntlnc,index_of_process_conv_trans)
              if(idtend>0) then
                dtend(:,:,idtend) = dtend(:,:,idtend) + gq0(:,:,ntlnc) - save_lnc
              endif
              idtend = dtidx(100+ntinc,index_of_process_conv_trans)
              if(idtend>0) then
                dtend(:,:,idtend) = dtend(:,:,idtend) + gq0(:,:,ntinc) - save_inc
              endif
            endif
          endif

        else
          do k=1,levs
            do i=1,im
              gq0(i,k,ntcw) = clw(i,k,1) + clw(i,k,2)
            enddo
          enddo
        endif   ! end if_ntiw

      else
        do k=1,levs
          do i=1,im
            clw(i,k,1) = clw(i,k,1) + clw(i,k,2)
          enddo
        enddo
      endif   ! end if_ntcw

    end subroutine GFS_suite_interstitial_4_run

  end module GFS_suite_interstitial_4